library(Seurat)
## Loading required package: ggplot2
## Loading required package: cowplot
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
## Loading required package: Matrix
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
seurobj <- readRDS('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/10x-180831-noreg')
#aligned <- readRDS('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/10x-180831-aligned')
#aligned_T1_old <- readRDS('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/10x-180504-T1-aligned')

Quality control plots

For now the data was only filtered on 0.1 percent.mito. The subcluster from T1 with low UMI and gene counts was removed.

VlnPlot(seurobj, c("nGene", "percent.mito", "nUMI"), group.by='timepoint', nCol = 1, point.size.use=-1, size.x.use = 10)

GenePlot(seurobj, 'nUMI', 'nGene', cex.use = 0.5)

PCElbowPlot(seurobj, num.pc=50) #TSNE+clustering run on 20 PC's.

Data visualization

Interesting to see: T4 and T5 contain a lot more variation than T1, T2 and T3, and PC2 seems to split T4 and T5. Could the split in PC2 describe the cells developing into white or brown?

PCAPlot(seurobj, group.by='timepoint', pt.size=0.1)

A few clusters in the data have much higher expression of ‘ADIPOQ’, ‘SCD’, ‘RBP4’, ‘G0S2’, ‘PLIN4’, ‘FABP5’. This is captured by PC2.

FeaturePlot(seurobj, reduction.use='pca', features.plot=c('ADIPOQ', 'SCD', 'RBP4', 'G0S2', 'PLIN4', 'FABP5'), pt.size=1, cols.use=c('gray', 'blue'), no.legend=F, nCol=2)

FeaturePlot(seurobj, reduction.use='pca', features.plot=c('PLA2G2A', 'MT1X', 'APOD', 'DPT', 'PTGDS', 'IGF2'), pt.size=1, cols.use=c('gray', 'blue'), no.legend=F, nCol=2)

PLA2G2A: http://www.jlr.org/content/early/2017/06/29/jlr.M076141 “…suggesting that PLA2G2A activates mitochondrial uncoupling in brown adipose tissue.”

PDGFRα/PDGFRβ signaling balance modulates progenitor cell differentiation into white and beige adipocytes. Based on PDGFRα or PDGFRβ deletion and ectopic expression experiments, we conclude that the PDGFRα/PDGFRβ signaling balance determines progenitor commitment to beige (PDGFRα) or white (PDGFRβ) adipogenesis. Our study suggests that adipocyte lineage specification and metabolism can be modulated through PDGFR signaling. http://dev.biologists.org/content/145/1/dev155861.long

FeaturePlot(seurobj, reduction.use='pca', features.plot=c('PDGFRA', 'PDGFRB'), pt.size=1, cols.use=c('gray', 'blue'), no.legend=F, nCol=2)

TSNE

TSNEPlot(seurobj, group.by='timepoint', pt.size=0.1)

TSNEPlot(seurobj, group.by='Phase', pt.size=0.1)

Cluster 11 = mixture cluster.

TSNEPlot(seurobj, group.by='res.0.5', pt.size=0.1, do.label=T)

VlnPlot(seurobj, group.by='res.0.5', features.plot=c('MALAT1', 'NEAT1'), point.size.use=-1)

Metadata feature plots

FeaturePlot(seurobj, reduction.use='tsne', features.plot = 'nUMI', cols.use=c('grey', 'blue'), no.legend=F)

FeaturePlot(seurobj, features.plot = 'percent.mito', cols.use=c('grey', 'blue'), no.legend = F)

FeaturePlot(seurobj, reduction.use='tsne', features.plot = 'nGene', cols.use=c('grey', 'blue'), no.legend = F)

Markergene expression

FeaturePlot(seurobj, reduction.use='tsne', features.plot = 'EBF2', cols.use=c('grey', 'blue'), no.legend = F)

FeaturePlot(seurobj, reduction.use='tsne', features.plot = 'TM4SF1', cols.use=c('grey', 'blue'), no.legend = F)

FeaturePlot(seurobj, reduction.use='tsne', features.plot = 'LY6K', cols.use=c('grey', 'blue'), no.legend = F)

FeaturePlot(seurobj, reduction.use='tsne', features.plot = 'PDGFRA', cols.use=c('grey', 'blue'), no.legend = F)

Markers literature

Marker genes for mature brown/beige compared to white mentioned by Seale 2016: UCP1, DIO2, CIDEA, PPARGC1A, PPARA, COX7A1, COX8B, PRDM16, EBF2. \

VlnPlot(seurobj, features.plot=c('UCP1', 'DIO2', 'CIDEA', 'PPARGC1A', 'PPARA', 'COX7A1', 'PRDM16', 'EBF2'), group.by='timepoint', point.size.use = -1, nCol=2)

Based on PDGFRα or PDGFRβ deletion and ectopic expression experiments, we conclude that the PDGFRα/PDGFRβ signaling balance determines progenitor commitment to beige (PDGFRα) or white (PDGFRβ) adipogenesis. Our study suggests that adipocyte lineage specification and metabolism can be modulated through PDGFR signaling. http://dev.biologists.org/content/145/1/dev155861.long

FeaturePlot(seurobj, reduction.use='pca', features.plot=c('PDGFRA', 'PDGFRB'), pt.size=1, cols.use=c('gray', 'blue'), no.legend=F, nCol=2)

FeaturePlot(seurobj, reduction.use='tsne', features.plot=c('PDGFRA', 'PDGFRB'), pt.size=1, cols.use=c('gray', 'blue'), no.legend=F, nCol=2)

GenePlot(SetAllIdent(seurobj, id='timepoint'), gene1='PDGFRA', gene2='PDGFRB', cex.use=0.5)

Stressed cells?

stress_genes <- read.table('/raid5/projects/timshel/sc-arc_lira/src/data-genelists/171219-van_den_Brink2017-genes_affected_by_dissociation.csv', header=T) %>% pull(1)
de_genes <- read.table('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/markergenes/180831/markers_res.0.7_negbinom', header=T)
de_genes <- de_genes[de_genes$p_val_adj < 0.05,]

Cluster 14 is the mixture cluster. Check MALAT1 and NEAT1 to be sure:

VlnPlot(seurobj, group.by='res.0.7', features.plot=c('MALAT1', 'NEAT1'), point.size.use =-1)

Are any of the stress genes DE genes for the mixture cluster?

de_genes_pos <- de_genes[de_genes$avg_logFC > 0, ]
de_genes_neg <- de_genes[de_genes$avg_logFC < 0, ]

intersect_pos <- list()
intersect_neg <- list()
n_pos <- list()
n_neg <- list()
perc_pos <- list()
perc_neg <- list()

for (i in unique(de_genes$cluster)){
  c <- paste('cluster', i)

  genes_pos <- de_genes_pos[de_genes_pos$cluster == i, 'gene']
  genes_neg <- de_genes_neg[de_genes_neg$cluster == i, 'gene']
  
  intersect_pos[[c]] <- length(intersect(toupper(stress_genes), genes_pos))
  intersect_neg[[c]] <- length(intersect(toupper(stress_genes), genes_neg))
  n_pos[[c]] <- length(genes_pos)
  n_neg[[c]] <- length(genes_neg)
  perc_pos[[c]] <- (intersect_pos[[c]] / n_pos[[c]]) * 100
  perc_neg[[c]] <- (intersect_neg[[c]] / n_neg[[c]]) * 100
}

data.frame(
  shared.pos.genes=unlist(intersect_pos),
  shared.neg.genes=unlist(intersect_neg),
  n_pos=unlist(n_pos),
  n_neg=unlist(n_neg),
  perc_pos=unlist(perc_pos),
  perc_neg=unlist(perc_neg)
)
mixture_genes_pos <- de_genes_pos[de_genes_pos$cluster == 13, ]
mixture_genes_neg <- de_genes_neg[de_genes_neg$cluster == 13, ]
intersect(toupper(stress_genes), mixture_genes_pos$gene)
## [1] "ACTG1"    "HSP90AB1"
VlnPlot(seurobj, features.plot=intersect(toupper(stress_genes), mixture_genes_pos$gene), point.size.use=-1, group.by='res.0.7', nCol=2)

VlnPlot(seurobj, features.plot='TXN', point.size.use=-1, group.by='res.0.7', nCol=2)

Gene correlations

matrix <- as.matrix(seurobj@data)
get_gene_correlations <- function(gene){
  gene <- as.numeric(matrix[gene,])
  correlations <- apply(matrix,1,function(x){cor(gene,x)})
  correlations <- as.data.frame(correlations[order(-correlations)])
  correlations['gene'] <- rownames(correlations)
  rownames(correlations) <- NULL
  names(correlations) <- c('cor', 'gene')
  return(correlations[order(-correlations$cor),])
}

PDGFRA. Positive correlations. (better to check per timepoint what the correlations are…)

pdgfra_correlations <- get_gene_correlations('PDGFRA')
## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero
pdgfra_correlations

PDGFRA negative correlations.

pdgfra_correlations[order(pdgfra_correlations$cor),]

Interesting: FABP5 as top negative correlated gene. Some of the other known adipogenic genes are also in the top (CIDEC, PLIN). Seale et al. describe PDGFRA as a general marker for preadipocytes and ADIPOQ (also top neg gene) as general marker for mature adipocytes. Then it makes sense that it is negatively correlated with the mature adipocyte genes.
Another study (2018) says:
Based on PDGFRα or PDGFRβ deletion and ectopic expression experiments, we conclude that the PDGFRα/PDGFRβ signaling balance determines progenitor commitment to beige (PDGFRα) or white (PDGFRβ) adipogenesis. Our study suggests that adipocyte lineage specification and metabolism can be modulated through PDGFR signaling. http://dev.biologists.org/content/145/1/dev155861.long
And:
In perigonadal (visceral) fat of male mice, beige adipocytes develop from a population of precursors that also differentiate into white adipocytes. These bipotent precursors express platelet-derived growth factor receptor-a. When the precursors proliferate, they lose PDGFRa expression and differentiate into UCP1+ adipocytes. Conversely, a high-fat diet stimulates the differentiation of PDGFRa+ cells into white adipocytes. This result is consistent with the finding that most or all white adipocytes descend from PDGFRa-expressing cells. Cell culture analyses have shown that PDGFRa-expressing cells can give rise to both UCP1- and UCP1+ (beige) adipocytes. 1 (do brown cells also express PDGFRa?)
Does that mean that PDGFRA is a general preadipocyte marker. If preadipocytes start to express PDGFRB they will develop into white and PDGFRA exprssion drops. Otherwise they will develop into beige and PDGFRA expression will also drop. So PDGFRA just serves to commit preadipocytes to beige and then it will drop? This seems to be in line with what we see in our data:

VlnPlot(seurobj, group.by='timepoint', features.plot=c('PDGFRA', 'PDGFRB'), point.size.use=-1, x.lab.rot=T)

EBF2 positive correlations.

ebf2 <- get_gene_correlations('EBF2')
## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero

## Warning in cor(gene, x): the standard deviation is zero
ebf2

Ebf2 negative correlations

ebf2[order(ebf2$cor),]